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ABSTRACT 

Faraday rotation of linearly polarized radio signals provides a very sensitive probe of fluctuations in 
the interstellar magnetic field and ionized gas density resulting from magnetohydrodynamic (MHD) 
turbulence. We used a set of statistical tools to analyze images of the spatial gradient of linearly 
polarized radio emission (|VP|) from the ISM for both observational data from a test image of the 
Southern Galactic Plane Survey (SGPS) and isothermal simulations of MHD turbulence. We compared 
the observational data with results of synthetic observations obtained with the simulations of 3D 
turbulence. Visually, in both data sets, a complex network of filamentary structures is seen. Our 
analysis shows that the filaments in the gradient can be produced by shocks as well as random 
fluctuations characterizing the non-differentiable field of MHD turbulence. The latter dominates for 
subsonic turbulence, while the former dominates for supersonic turbulence. We show that these 
two regimes exhibit different distributions as well as different morphologies in the maps of |VP|. 
Particularly, filaments produced by shocks show a characteristic "double-jump" profile at the sites of 
shock fronts resulting from delta function like increases in the density and/or magnetic field. Filaments 
produced by subsonic turbulence show lower values of |VP|, have a single jump profile, and are 
more dominated by fluctuations in the magnetic field. In order to quantitatively characterize these 
differences we use the topology tool known as a genus curve as well as the moments of the image 
distribution, i.e. the mean, variance, skewness and kurtosis. We find that higher values for the 
moments correspond to cases of |VP| with larger Mach numbers, but the strength of the dependency 
is connected to the telescope angular resolution: higher resolution yields stronger differences between 
the sonic Mach number regimes. In regards to the topology, the supersonic filaments observed in 
I VP I have a positive genus shift, which indicates a "swisscheese" like topology, while the subsonic 
cases show a negative genus, indicating a "clump" like topology. Transonic simulations have a more 
neutral topology. In the case of the genus, the dependency on the telescope resolution is not as strong, 
making this tool particularly useful for observational studies. Based on the analysis of the genus and 
the higher order moments, the SGPS test region data has a distribution and morphology that matches 
subsonic to transsonic type turbulence, which independently confirms what is now expected for the 
WiM. Our initial statistical application to gradients of linear polarized radio data shows this type 
of data to be very promising for studies of turbulence in the ISM in the way of obtaining the Mach 
numbers. 

Subject headings: ISM: structure — MHD — turbulence. Linear Polarization, Faraday Rotation 



1. INTRODUCTION 

The interstellar medium (ISM) of our Milky Way is 
host to a variety of physical mechanisms that regulate 
and govern the structure and evolution of the Galaxy. 
The current understanding of the ISM is that it is a 
multi-phase environment composed of a tenuous plasma, 
consisting of gas and dust, which is both magnetized and 
highly turbulent (Ferriere 2001, McKee & Ostriker 2007). 
In particular, the awareness of turbulence as a dominant 
physical process in the ISM has only happened in the 
last decade (Elmegreen & Scalo 2004). Turbulence plays 
a critical role in the areas of star formation, magnetic 
reconnection, magnetic field amplification, nearly every 
transport process, cosmic ray acceleration, magnetic dy- 
namo, and the physics in the intercluster medium, to 
name just a few (see Lazarian & Vishniac 1999, Vishniac 
& Cho 2001, Elmegreen & Scalo 2004, Lazarian 2006, 
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Ballesteros-Paredes et al. 2007, McKee & Ostriker 2007 
and references therein) . Additionally, turbulence has the 
unique ability to transfer energy over scales ranging from 
kiloparsecs down to the proton gyroradius. This is criti- 
cal for the ISM, as it explains how energy is distributed 
from large to small spacial scales in the Galaxy. 

Despite the now obvious importance of magnetized tur- 
bulence, the situation of understanding ISM physics is no 
less complicated. In spite of the big recent advances in 
understanding of incompressible and compressible MHD 
turbulence (see Goldreich & Sridhar 1995, Cho & Vish- 
niac 2000, Maron & Goldreich 2001, Cho, Lazarian & 
Vishniac 2002, 2003, Cho & Lazarian 2003, Kowal & 
Lazarian 2010, Beresnyak & Lazarian 2010, Beresnyak 
2011) the ISM presents a complex environment with mul- 
tiple energy injection sources, different phases and vari- 
ous instabilities acting at different scales. The properties 
of this turbulence affect key astrophysical processes and 
obtaining properties of ISM turbulence from observations 
opens ways of gauging numerical simulations and testing 
theory. In light of these complexities, the most fruitful 
way of studying astrophysical MHD turbulence and the 
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Fig. 1. — Top: Observational P (left) and |VP| (right) from tiie SGPS data used in this study. Bottom: MHD simulation model number 
1 in Table [Tl with P (left) and |VP| (right). The simulated map of P has had the mean subtracted from maps of Q and U, similar to the 
observations. 



processes it affects is to use a synergetic approach which 
combines the knowledge of theoretical predictions, nu- 
merical studies, and observational efforts. 

Observationally there are several ways of studying 
MHD turbulence. Many of these techniques hinge on 
density fluctuations in ionized or neutral media (see 
Spangler & Gwinn 1990, Armstrong et al. 1995, Padoan 
et al. 2003, Falgarone et al. 2005, Chepurnov & Lazar- 
ian 2010) and are aimed at finding the density power 
spectrum of turbulence. More recently, ways to find 
the magnetization and the Mach number of turbulence 
have been explored (Kowal, Lazarian & Beresnyak 2007, 
Burkhart et al. 2009, 2010, Esquivel & Lazarian 2010, 
Tofflemire, Burkhart & Lazarian 2011, Esquivel & Lazar- 
ian 2011). While column density data are arguable the 
most common type, they do not contain the full 3D pic- 
ture and are only passive tracers of the turbulence veloc- 
ity field. Various ways to study the interstellar velocity 
field have been explored. The tested and theoretically 
motivated ways to study turbulent velocities in super- 
sonic interstellar turbulence are based on the use of the 
Velocity Channel Analysis (VGA) and Velocity Coordi- 
nate Spectrum techniques (Lazarian & Pogosyan 2000, 
20004, 2006, 2008). These techniques have been used 
both for atomic HI and molecular data (see Stanimirovic 
& Lazarian 2001, Padoan et al. 2006, 2009, Chepurnov 
et al. 2010) to find the spectra of the turbulent velocity 
fields. 

Naturally, studies of turbulence and magnetic fields are 
of great importance and synergetic value for the ISM. In 



this case, a number of techniques have been explored, 
e.g. structure functions of the polarization vectors aris- 
ing from dust polarized emission (see Falceta-Gongalves 
et al. 2009, Houde et al. 2011). The quantitative study 
of synchrotron intensity fluctuations can be traced to 
works by Getmansev (1958), while fluctuations of syn- 
chrotron polarizatior0 were used, for instance, to evalu- 
ate the spectra of magnetic turbulence in Hydra cluster 
(Enfilin & Vogt 2006, Enfihn et al. 2010) 

More recently, several authors have discussed the 
prospects of the use of radio polarization maps to study 
turbulence (see Haverkorn & Heitsch 2004, Fletcher & 
Shukurov 2006, 2007, Gaensler et al. 2011). Faraday ro- 
tation maps of linearly polarized radio signals are espe- 
cially promising as they provide very sensitive probes of 
fluctuations in magnetic field and ionized gas density (see 
Gray et al. 1998, Gaensler et al. 2001, and Landecker et 
al. 2010). The Faraday rotation can be calculated as: 



RM = K 



.(/)B(l)dl 



(1) 



With units of rad m^^, K = 0.81 rad m^^ pc^^ cm'^ 
IJ,G~^ and B, I and rie are the magnetic field strength 
in fiG, the distance along the LOS in parsecs, and the 
electron density in crn~"^ along the LOS, respectively. 

^ For a theoretical description of synchrotron fluctuations for 
arbitrary index of cosmic rays and realistic models of anisotropic 
turbulence see Lazarian &; Pogosyan (2011). 
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Although many objects seen in the polariza- 
tion/Faraday maps can be matched to objects seen in 
other wavelengths (such as supernova remnants), an ex- 
tended diffuse polarization emission network that is rich 
in structure is also present that can not be mirrored 
in other wavebands or in total intensity (Fletcher & 
Shukurov 2007) . The intensity variations seen in maps of 
Q, U and P are the result of small-scale angular structure 
in the Faraday rotation induced by foreground ionized 
gas and are thus an indirect representation of turbulent 
fluctuations in free electron density and magnetic field 
throughout the ISM. 

In this paper we use polarization gradients to study 
ISM turbulence. The use of gradients to highlight small 
rapid fluctuations seen in polarization maps was first dis- 
cussed by Gaensler et al. (2011). When the spatial gra- 
dient is applied to maps of vector P — {Q,U) a com- 
plex web of filamentary structures is revealed. These 
filaments (see right column of Figure [1] for an example) 
were interpreted by Gaensler et al. (2011) as rapid fiuc- 
tuations in Ue and B along the LOS due to turbulence. 
In this paper we will further explain the origins of the 
filamentary structures as they are related to turbulence 
and develop quantitative methods that can be applied to 
this data in order to obtain the Mach numbers of turbu- 
lence. Taking the gradient of rotation measure or linear 
polarization maps has its advantages and disadvantages. 
The primary advantage is that the spatial gradient of 
P satisfies the property that it has both translational 
and rotational invariance in the (Q,U) plane. Quanti- 
ties such as the polarization amplitude and polarization 
angle are not preserved under arbitrary translations and 
rotations, which can result from one or more of a smooth 
distribution of intervening polarized emission, a smooth 
uniform screen of foreground Faraday rotation, or the 
effects of missing large-scale structure in an interfero- 
metric data-set. Thus the magnitude of the gradient is 
the simplest quantity that is not significantly effected by 
missing large-scale structure or excess foreground emis- 
sion or Faraday rotation. Taking the gradient allows 
one to clearly see jumps and discontinuities, regardless 
of whether single-dish (total power) measurements are 
present in the data. In particular, this will highlight ar- 
eas where a sharp change in Ue or B occurs, which is 
most likely due to turbulent fluctuations or shock fronts 
in the ISM. However, one must also keep in mind the dis- 
advantages of using gradients, namely the fact that the 
gradient may enhance noise. 

In this work we explore the physical causes of the fil- 
aments by taking the gradient of polarization maps of 
isothermal MHD turbulence. We investigate the depen- 
dency of the sonic Mach number on the structures seen 
in |VP| and calculate measures of the Probability Distri- 
bution Functions (PDFs) as well as the a measure of the 
topology called genus. Genus has been extensively used 
in cosmology studies (Gott et al. 1986) and has been 
suggested for ISM studies in Lazarian (1999). Its use for 
synthetic column density maps and observations is dis- 
cussed in the literature (Lazarian, Pogosyan & Esquivel 
2002, Kowal, Lazarian & Beresnyak 2007, Kim & Park 
2007, Chepurnov et al. 2008). 

The PDFs and genus have been studied on MHD tur- 
bulence in the past and have both shown sensitivity 
to the sonic Mach number in density, column density 




Fig. 2. — Examples of LOS maps and their respective gradients 
relevant to this paper for subsonic turbulence (model 1). The first 
column shows column density, LOS magnetic field, and the rota- 
tion measure from top to bottom. The second column shows the 
gradients of column density, LOS magnetic field, and polarization 
vector. 

and position-postion-velocity data (Padoan et al. 1999, 
Kowal, Lazarian & Beresnyak 2007, Chepurnov et al. 
2008, Burkhart et al 2009, 2010, Tofflemire et al. 2011). 
We view these statistics as part of a set of tools which 
can be applied to polarization data in order to determine 
the sonic Mach number of the turbulence in the ionized 
ISM. We stress that the use of the whole set provides a 
synergetic quantity: obtaining the same result with dif- 
ferent techniques substantially increases how trustworthy 
the result is. 

The paper is organized as follows. In § [5] we further 
describe the data sets used, in particular the Southern 
Galactic Plane Survey and a set of ideal MHD simula- 
tions, and our calculation of the rotation measure, the 
linear polarization maps and their gradients. In § [3] we 
discuss the origins of the observed filaments as they re- 
late to the sonic Mach number. We describe different 
statistical measures of the sonic Mach number in § [U in 
particular the genus and PDFs. In we discuss our results 
followed by conclusions in § [5l 

2. DATA AND METHOD 

2.1. Gradient Technique and relation of\VP\ to \VRM\ 

An important observational quantity connected to in- 
terstellar density and magnetic field fluctuations is the 
Faraday effect. In the presence of magnetic fields and free 
electrons, bifringence of circularly polarized orthogonal 
modes occurs, giving these modes two different propaga- 
tion velocities. In the case of pure polarized background 
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emission propagation through a magnetoionized medium, 
the hnearly polarized radiation will emerge with its po- 
larization position angle rotated by the amount given 
in Equation [T] Thus the relation between the observed 
position angle, emitted position angle, and the rotation 
measure is: 

e - Go = RMX^ (2) 

which has units of radians. Here A is the wavelength in 
meters. 



VN 




Fig. 3. — Examples of LOS maps and their respeetive gradients 
relevant to this paper for supersonic turbulence (model 6). The 
first column shows column density, LOS magnetic field, and the 
rotation measure from top to bottom. The second column shows 
the gradients of column density, LOS magnetic field, and polariza- 
tion vector. 

Observational determination of Faraday rotation 
comes from measurements of the linear polarization vec- 
tor P = (Q,U) (which depends on Stokes U and Q as 

\P\ — \J ^ U^) as a function of A^. To avoid confu- 
sion between vector and scalar P, we use bold notation 
to denote the vector quantity of the linear polarization 
map. 

We define the gradient of the polarization vector as: 

(3) 

We note that in the case of vector P we have: 



From this equation, one can derive a relationship be- 
tween |VP| and |Vi?M| for Faraday-thin polarized emis- 
sion as: 

\VRM\ = |VP|/2iA2|P| (5) 

Equation 5 only holds for data in which the entire sig- 
nal is measured (i.e. singic-dish data included) and for 
which the background is uniform. When |P| = 1 and 
A = 1 (which are the assumptions we use for the sim- 
ulations), one finds a trivial relation between |VP| and 
\VRM\ as \VRM\ = |VP|/2i. However, this relation can 
only be used to calculate |VP| or |Vi?Af| in the simula- 
tions, since the assumption of |P| = 1 is almost always 
to simplistic for the observations because the data are 
missing single-dish information and/or the background 
I VP I is not zero. 

We calculate maps of RM, P and their gradients from 
density and LOS magnetic field that is perpendicular and 
parallel to the mean magnetic field in the simulations. 
We calculate the RM as per Equation [T] at every point 
and then take its spatial gradient, that is, we compute the 
gradient vector at every pixel of the image using neighbor 
pixels. We can calculate |P| by calculating the stokes 
vectors as Q = cos(20) U = sm{29). These expressions 
come from applying Equation 2 with assumed values for 
A and 9q. 

We show a subsonic and supersonic case of density (n) , 
Vn, LOS magnetic field (LOS B), VS, RM and |VP| in 
Figure [2] and Figure [3l respectively. A comparison of the 
SGPS test data and a subsonic case is given in Figure 
[TJ Inspection of maps of |VP| reveals that filaments are 
created in both cases and that there is some correlation 
between gradients of column density, magnetic field, and 
|VP|. We will discuss these further in Section[3l 

2.2. Southern Galactic Plane Survey 

We use a subsection of radio continuum images of an 
18-square-degree patch of the Galactic plane, observed 
with the Australia Telescope Compact Array (ATCA, see 
McClure-Griffiths et al. 2001 and Gaensler et al. 2001 
for more details). We examine the 1.4 GHz frequency 
data averaged over adjoining frequency channels with si- 
multaneously recorded Stokes /, Stokes Q, and Stokes U 
as part of the SGPS test region (Gaensler et al. 2001). 
This field consists of 190 mosaicked pointings of the Aus- 
tralia Telescope Compact Array (ATCA) and covers the 
range 325.5 < / < 332.5, -0.5 <b < 3.5. Complicated 
extended structure is seen in linear polarization through- 
out the test region, almost all of which has no correlation 
with total intensity. We select a 512 x 512 pixel subregion 
from this data to match the resolution of the simulations 
used in our study and display it in Figure [T] in the top 
row. The SGPS region we select begins at coordinate 
1=332.3373, b= -0.3138 and is not overly contaminated 
by bad pixels and contains significant emission. 

2.3. Simulations 

We generate a database of 3D numerical simulations of 
isothermal compressible (MHD) turbulence by using the 
MHD code of Cho & Lazarian (2003) and varying the 
input values for the sonic and Alfvenic Mach number. 
The sonic Mach numbers are defined as Ms = (|v|/Cs), 
where is v is the local velocity, Cg is the sound speed, 
and the averaging is done over the whole box. Similarly, 
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TABLE 1 

Description of the simulations 



MHD, 512'^ 



iVioaei 


Pgas 


i>oxt 


JVls 




Description 


1 


2.00 


1.00 


0.5 


0.7 


subsonic &z sub-Alfvcnic 


2 


0.70 


1.00 


1.0 


0.7 


transsonic & sub-Alfvcnic 


3 


0.10 


1.00 


2.0 


0.7 


transsonic & sub-Alfvcnic 


4 


0.05 


1.00 


3.0 


0.7 


supersonic & sub-Alfvenic 


5 


0.025 


1.00 


4.4 


0.7 


supersonic & sub-Alfvenic 


6 


0.0077 


1.00 


8.0 


0.7 


supersonic & sub-Alfvenic 


7 


0.0049 


1.00 


10 


0.7 


supersonic & sub-Alfvenic 


8 


2.00 


0.10 


0.5 


2.0 


subsonic & super-Alfvenic 


9 


0.70 


0.10 


1.0 


2.0 


transsonic & super-Alfvenic 


10 


0.10 


0.10 


2.0 


2.0 


transsonic & super-Alfvenic 


11 


0.05 


0.10 


3.0 


2.0 


supersonic & super-Alfvenic 


12 


0.025 


0.10 


4.4 


2.0 


supersonic & super-Alfvenic 


13 


0.0077 


0.10 


8.0 


2.0 


supersonic & super-Alfvenic 


14 


0.0049 


0.10 


10 


2.0 


supersonic & super-Alfvenic 



the Alfvenic Mach number is A4a = (IvI/wa), where 
VA — |B|/y/p is the Alfvenic velocity, B is magnetic field 
and p is density. We briefly outline the major points of 
the numerical setup (for more details see Cho & Lazarian 
(2003). 

The code is a second-order- accurate hybrid essentially 
nonoscillatory (ENO) scheme (Cho & Lazarian 2003) 
which solves the ideal MHD equations in a periodic box: 



§ + V-(pv)=0, (6) 



dpv 



V • 



pvv 



-BB 



= f, (7) 



— -Vx(vxB)=0, (8) 

with zero-divergence condition V • B = 0, and an isother- 
mal equation of state p = C^p, where p is the gas pres- 
sure. On the right-hand side, the source term f is a 
random large-scale solenoidal driving forceQ. 

The magnetic field consists of the uniform background 
field and a fluctuating field: B — Bcxt + b. Initially 
b = 0. 

We scale the simulations to physical units, adopting 
typical parameters for warm ionized gas. We assume a 
pixel size of 0.15 parsecs and density of 0.1 cm^"^. The 
simulations are assumed to be fully ionized and we do 
not include the effects of partial ionization. 

To make the maps of |VP| we first calculate the LOS 
rotation measure at each pixel then we take the take 
the gradient of this rotation measure map and convert 
it to |VP| via Equation [5l An equally valid way is to 
calculate the rotation measure at each pixel, then shine 
a polarized signal through the cube with our assumed 
background values of Q = 1, U = 0. Then one can 
calculate the emergent values of Q and U by applying 
the simulated rotation measure map and then calculate 
the resulting |VP|. Both methods will produce identical 
maps of |VP|. Additional smoothing of the maps of Q 
and U using a Gaussian kernel can also be performed to 
mimic the telescope beam. 

4 f = pdv/dt 



3. THE ORIGIN OF FILAMENTARY STRUCTURES IN 
POLARIZATION AND ROTATION MEASURE 
GRADIENTS 

A filament traced by |VP| or |Vi?Af| will form as a 
result of a localized change in either density or magnetic 
field as a function of position on the sky as per Equation 
[TJ There are many physical processes that can result 
in sharp changes in these quantities in the ISM, includ- 
ing gravitational collapse and outflows. However, high 
Reynolds fluids in the ISM are expected to be turbulent 
(see Cho, Lazarian & Vishniac 2003, Elmegreen & Scalo 
2004, Lazarian et al. 2009) and a more ubiquitous pro- 
cess responsible for fluctuations in or B is due to MHD 
turbulence, precisely because it expected everywhere in 
the ISlVfl, although the type of turbulent environment 
can vary (e.g. the compressibility, magnetization, equa- 
tion of state, etc.). 

In the ISM, fluctuations in density and magnetic field 
will occur as a result of MHD turbulence, which will be 
visible in polarimetric maps. In the case of taking gra- 
dients of a turbulent field, one would expect to find fila- 
mentary structure created by shock fronts, jumps and 
discontinuities. Figure |4] shows a cartoon illustrating 
these three separate cases of a possible profile and its 
respective derivative. The cases are: 

1. A Holder continuous profilcH that is not differen- 
tiable at a given point (e.g. the absolute value 
function at the origin): Common for all types of 
MHD turbulence. 

2. A jump profile: Weak shocks, strong fluctuations 
or edges (e.g. a cloud in the foreground which sud- 
denly stops). 

3. A spike profile (e.g. delta function): Strong shock 
regime. 

In respect to case one, it is known that the turbulent 
velocity field in a Kolmogorov-type inertial range both in 
hydro and MHD are not differentiable, but only Holder 
continuous (Bernard et al. 1998, Eyink 2009). Another 
example is that of any fractal function that displays self- 
similarity but is not differentiable everywhere. This pro- 
file will naturally create discontinuities when one takes 
its derivative. Therefore, case one can be found in both 
subsonic and supersonic type turbulence. Case one type 
filaments can be seen in Figure [2] for N, B, and P in the 
right column. Case two creates a structure in the gra- 
dient by a shock jump or a large fluctuation in either 
rie or B. Here again, this type of enhancement in |VP| 
could be found in supersonic and subsonic type turbu- 
lence and is due either to large random spatial increases 
or decreases due to turbulent fluctuations along the LOS 
or weak shocks. We expect weak shock turbulence to 
show a lager amplitude in |VP| then the subsonic case. 
Case three is unique to supersonic turbulence in that it 
represents a very sharp spike in rig and/or B across a 
shock front. The difference between this case and what 

^ Because of the large injection scale of the ISM, the Reynolds 
numbers can typically reach 10^'' 

^ Holder continuous functions satisfy \u{xi,t) — u(x2,t)\ < 
C\xi — X2^ ■ When the exponent h=l, this satisfies the Lipschitz 
condition. In the case of turbulent fields h = \^ 



6 



BURKHART ET AL. 



might be seen in case two is that here we are deahng with 
interactions of strong shock fronts, which are known to 
create deha function hke distributions in density (Kim 
& Ryu 2005). In this case, the derivative of case three is 
has a distinctly different profile with respect to case one 
and two. Case three shows a 'double jump' profile across 
the shock front, which can be seen in Figure [3] in the top 
right panel for LOS density and the bottom right panel 
for |VP|. This morphological distinction can be used to 
determine if one is dealing with turbulence that is in a 
shock dominated regime (i.e. supersonic) and can pro- 
vide researchers with a promising new avenue of obtain- 
ing the sonic Mach numbers from polarimetric data. In 
the case of the Alfvenic Mach number, the morpholog- 
ical difference is less clear, however gradients will tend 
to align along the field lines in the case of strong field 
(sub- Alfvenic turbulence). 



Case 1 



dn/dx 




Case 2 



dn/dx 



Case 3 



dn/dx 



Fig. 4. — Cartoon example of three possible scenarios for en- 
hancements in a generic image "n" , where "n" could be P|, RM 
or p (density or EM). Case one (top row) shows an example of 
a Holder continuous function that is not diffcrentiable at the ori- 
gin (applicable to all turbulent fields). Case two (middle row) 
shows an example of a jump resulting from strong turbulent fluc- 
tuations along the LOS or weak shocks, case three (bottom row) 
shows a delta function profile resulting from interactions of strong 
shocks. In this case, the derivative gives a 'double jump' profile 
which produces morphology that is distinctly different from the 
previous cases. 

Also of interest is the question of which quantity is 
providing the dominate contribution to the structures in 
|VP| or \VRM\: Vn^, VB^os or both equally? Espe- 
cially in the case of compressible turbulence, the mag- 
netic energy is correlated with density, namely, denser 
regions contain stronger magnetic fields, which is due 
to the compressibility of the gas (Burkhart et al. 2009). 
This causes the magnetic field to follow the fiow of plasma 
if the magnetic tension is negligible. The compressed 



regions are dense enough to distort the magnetic field 
lines, enhance the magnetic field intensity, and effectively 
trap the magnetic energy due to the frozen-in condition. 
Thus, for the supersonic cases, the intensity of the struc- 
tures seen in VP are more pronounced then in the sub- 
sonic case, which is observed when comparing Figures [2] 
and [21 

However, in the case of subsonic turbulence, there are 
no compressive motions. In this case, random fluctua- 
tions in density and magnetic field will create structures 
in |VP| and \V RM\. 

Due to these effects, we might expect different trends 
in the correlation of supersonic and subsonic |VP| with 
ViV or VEM (the gradient of the emission measure) and 
V-B. We test this by plotting the pixel by pixel corre- 
lation coefficient of |VP| with the gradients of EM, N, 
and LOS B in Figure [5] In the case of subsonic turbu- 
lence, I VP I better traces out the fiuctuations in \7B (blue 
line), while the supersonic cases are dominated by den- 
sity fluctuations. This is because density enhancements 
are dominate due to shock fronts in the case of super- 
sonic turbulence, while in subsonic turbulence density is 
highly incompressible. In this case, the magnetic field 
will dominate the topology of the rotation measure and 
|VP|. This behavior is analogous to velocity in neutral 
hydrogen radio position-position- velocity cubes of turbu- 
lence, where density dominates the power spectrum for 
the case of supersonic turbulence and velocity dominates 
the spectrum for subsonic turbulence (see Lazarian & 
Pogosyan 2006, Burkhart et al. 2011a). This difference 
in correlation provides yet another way of gauging the 
Mach numbers if one has both the emission measure and 
the linear polarization map. Correlated spatial gradients 
between the two should indicate regions of shocks. 

In the next section we will explore the utility of gra- 
dients of polarimetric data for the determination of the 
Mach numbers by investigating two different statistical 
measures of looking at the distribution and topology of 
the I VP I maps: PDF moments and genus function. 

4. STATISTICAL DETERMINATION OF THE SONIC MACH 

NUMBER 

The previous section provided some theoretical discus- 
sion for why we expect |VP| data to be useful for de- 
termining the sonic Mach number. In this section we 
will attempt to statistically quantify the differences seen 
in both the morphology and the distribution of maps of 
|VP|. We again note our assumption for the simulations 
of |P| = 1 thus giving a trivial scaling relationship be- 
tween |VP| and |Vi?M| as: |Vi?M| = |VP|/2. We also 
provide an observational comparison for both statistics 
with the SGPS test region shown in Figure [TJ 

4.1. Moments 

A probability distribution function (PDF) is the func- 
tion describing the frequency of occurrence of values in 
the distribution of intensities. PDFs and their quanti- 
tative descriptors have been used to study turbulence 
in a variety of astrophysical context including diffuse 
ISM turbulence (Berkhuijsen & Fletcher 2011), turbu- 
lence characterization (Federrath et al. 2010, Esquivel et 
al. 2010, Audit & Hennebelle 2010) solar wind (Burlaga 
et al. 2007) and molecular ISM (Padoan et al. 1999). 
Several authors have discussed the use of PDFs in de- 
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Fig. 5. — Correlation coefficient between|VP| and the LOS S7EM 
(black diamonds), VAf(red asterisk) and \7B (blue plus sign). The 
left panel is super-Alfvenic and the right panel is sub-Alfvenic. In 
the case of strong shocks, strong correlation is observed between 
S7EM and VA^ due to enhanced density. Subsonic models show 
strong correlations with V-B. 

termining the Mach numbers of ISM turbulence, in par- 
ticular the sonic Mach number (Vazquez-Semadeni 1994, 
Padoan et al. 1999, Kowal, Lazarian & Beresnyak 2007, 
Burkhart et al. 2010, Price, Federrath, & Brunt 2011). 
However, this technique is almost always used in the con- 
text of the column density. To our knowledge, no one has 
applied this technique to the rotation measure, polariza- 
tion maps, or their gradients. 

One method of describing PDFs is by using statisti- 
cal moments to characterize the mean and variance and 
departures from Gaussianity. 

The first and second order statistical moments (mean 
and variance) used here are defined as follows: /i^ = 

■h E^i and = ^ Xl^i {^i - ?)^ respectively. 

The 3rd and 4th order moments, Skewness and kurtosis 
respectively, are defined as: 




Where N is the total number of elements and is the 
distribution of intensities. 

Past works have focused on the relationship between 
the moments and the density or column density. As the 
Mach number increases, so does the mean value of den- 
sity, as shocks increase density (the sonic Mach num- 
ber goes as Ms ~ P^^^)- The variance, skewness, and 
kurtosis are less obvious quantities and their relation to 
the sonic Mach number has been derived numerically 
(Padoan et al. 1999, Kowal, Lazarian & Beresnyak 2007, 
Burkhart et al. 2010 Price, Federrath, & Brunt 2011). 
Variance has some disadvantages to skewness and kur- 
tosis. One is that variance is scale dependent, and val- 
ues will change between different data set normalizations. 
This makes direct comparison between simulations and 
observations difficult. On the other hand, the higher 
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order moments (skewness and kurtosis) describe devia- 
tions from Gaussianity and are unit-less numbers, and 
therefore are scale free. They also are shown to increase 
more linearly with the sonic Mach number in the case 
of column density (Kowal, Lazarian & Beresnyak 2007, 
Burkhart et al. 2010). However, the increase is not so 
pronounced in the case of subsonic turbulence, making 
variance a better indicator in this regime. 
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Fig. 6.— PDFs of |VP| normalized as f(x) = (|VP|- < |VP| > 
)/o-(|VP|) where (7(|VP|) is the standard deviation of |VP|. The 
supersonic case (represented by -|- signs) is highly skewed and kur- 
totic compared with the subsonic case (represented by * signs). 

While the rotation measure is obviously related to the 
density, the relationship between the moments of |VP| 
and the Mach number has never been studied. The po- 
larization gradient may have an advantage over column 
density/dispersion measure for investigating the Mach 
numbers of the WIM, as shocks will cause both an in- 
crease in density and magnetic field due to correlations 
between density and field strength in supersonic type tur- 
bulence (see Burkhart et al. 2009). Furthermore, the 
gradient will highlight interesting features in the obser- 
vational data which might be buried by observational 
effects, such as a DC offset, that otherwise would not be 
seen in the maps of column density or linear polarization 
(see Figures [5] and [3]) . 

Figure |6] shows an example PDF of a subsonic and su- 
personic simulation (model 1 and 6 ). Here it is clear that 
the general relations that are true for the column density 
will apply here: as the Mach number increases the distri- 
bution becomes more kurtotic and more skewed. Figure[7] 
shows the moments of |VP| vs. the sonic Mach number 
for our simulations with LOS taken along the mean mag- 
netic field. Error bars are created by taking the standard 
deviation between different time snapshots of the data. 
We see that all four moments tend to increase with in- 
creasing sonic Mach number. This is due to the shock 
fronts creating more discontinuities and sharper gradi- 
ents. This not only increases the mean value and the 
variations in |VP|, but also creates very peaked distri- 
butions and distributions with tails skewed towards the 
right. 

In terms of the Alfvenic Mach number, the main con- 
tributor is the LOS direction chosen relative to the or- 
dered field. When we look parallel to the mean field 
as shown in in Figure [Tj the sub-Alfvenic cases show 
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higher mean and variance while the super- Alfvenic cases 
show generally higher skewness and kurtosis. We plot 
the mean and variance of the polarization angle in Fig- 
ure [HI which confirms the trend seen in the gradients. 
In the case of the LOS parallel to the mean field, 
the mean value of the polarization angle has a strong 
Alfvenic dependency with no sonic Mach number de- 
pendency. The variance shows strong dependency on 
both sonic and Alfvenic Mach numbers. However, these 
trends are unique for sight-lines parallel to the mean field 
LOS. When looking perpendicular to the LOS the ef- 
fects change. For example, the moments of |VP| show 
stronger skewness and kurtosis in the case of the sub- 
Alfvenic models. Thus to gauge the Alfvenic Mach num- 
ber using the distribution, the ordered LOS magnetic 
field is required. 

An additional effect that must be considered is the 
issue of the telescope resolution. For example. Figure |9] 
shows model number 6 without smoothing (left) and with 
Gaussian smoothing with FWHM=6 pixels (right). It is 
clear that smoothing changes the distribution of maps of 
|VP|. We plot the moments vs. smoothing FWHM for 
four models in Figure [101 As the resolution of maps of 
Q and U decrease, so do the moments. Thus one needs 
to take the smoothing of the data into account when 
comparing PDFs of |VP|. 

We investigate the moments of regions of the SGPS 
data with the most emission removing the bad pixels 
from our investigation. We find that the SGPS data has 
average skewness of 0.825 and kurtosis of 0.928 and mean 
value and variance of 0.004 and 3.38a:10-^ respectively. 
Comparing with our numerical set up, this most closely 
matches subsonic to transonic values of the moments, 
with telescope smoothing taken into account. 

4.1.1. Moment Maps 

From section 4.1 it is clear that, for a given map of 
|VP|, as the sonic Mach number increases the skewness 
and kurtosis also increase. However, these were globally 
averaged values of both the moments (i.e. the PDF of 
the entire image) and the sonic Mach number. If we want 
to look at smaller scale variations, we can calculate the 
PDFs of smaller portions of the image. Using a moving 
kernel (box, circle etc.), we can create a "moment map" 
which is essentially a smoothed map that calculates the 
moments at every point with a given box size. Because 
we are dealing with gradient quantities, we expect these 
maps to particularly trace regions where shocks are inter- 
acting and along shock fronts, thus tracing areas where 
the sonic Mach number is changing. 

We make skewness and kurtosis moment maps of the 
|VP| images similar to the method of Burkhart et al. 
(2010). This results in a smoothed moment image, which 
can be compared to the actual image of the LOS local 
sonic Mach number (LOS average sonic Mach number 
for each pixel) with the same smoothing kernel as the 
|VP| moment map. In this section, we investigate the 
relationship between the LOS sonic Mach number and 
the values of skewness and kurtosis on a pixel-by-pixel 
basis. 

The key questions one must consider when making a 
moment map are what the kernel shape should be, what 
resolution is appropriate for good statistics, and how to 
compare the values of the higher order moments to the 



sonic Mach number? The last point is particularly im- 
portant. While Burkhart et al. (2010) took a linear 
fit between the sonic Mach number and the moments in 
the case of column density, we will test how effective 
the moment maps can distinguish between supersonic 
and subsonic regimes on a pixel-by-pixel bases. Thus 
we stress that we are not trying to calculate the exact 
sonic Mach number with this method, which is a topic 
for another work. We only seek to determine how well 
a moment map is able to pick out regimes of supersonic 
and subsonic. For our initial test, we choose to use a 
boxcar kernel for both the |VP| maps and for averaging 
the sonic Mach number, since a boxcar can handle edges 
in a square image most effectively. 




Fig. 7. — Moments of |VP| vs. Ms- Error bars are created by 
taking the standard deviation of the moments between different 
time snapshots of the well-developed turbulence. We show sub- 
Alfvcnic cases in red and supcr-Alfvenic cases in black. 

How well the moment map is able to determine if the 
gas is supersonic or subsonic depends on three parame- 
ters: the threshold level for skewness {jt, above which 
the gas is considered supersonic and below it is con- 
sidered subsonic), the threshold level for kurtosis (/3t, 
above which the gas is considered supersonic and below 
it is considered subsonic), and the kernel box size chosen. 
The values for jt and I3t represent the threshold value 
between supersonic and subsonic regimes. Any pixels be- 
low either 7t and is deemed subsonic and any pixel 
above both is deemed supersonic. Of course, an inter- 
mediate threshold value could also be chosen to probe 
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Fig. 8. — Mean (top) and variance (bottom) of the polarization 
angle (or RMA^) along the direction of the mean magnetic field vs. 
sonic Mach number. A strong Alfvenic dependency is observed. In 
the case of the mean value, the dependency is independent of the 
sonic Mach number. The trend is different for a different LOS. 

the transsonic regime, but we omit this here. Again we 
should stress that we are not employing this method to 
get an exact Mach number; we simply view it as a way 
of determining whether our |VP| image shows supersonic 
or subsonic characteristics. 

The best fitting parameter values of box size, jt and 
Pt are those that provide accurate translation between 
the moments and the LOS Mach number. In order to 
determine these values we use a genetic algorithm which 
searches possible combinations of these three parameters 
in order to determine the best fit. The genetic algorithm 
can provide more computationally efficient way of test- 
ing for fitness rather then iterating over every possible 
combination of parameter space. See the appendix for a 
detailed description of our algorithm. 

Although multiple combinations of parameters showed 
high confidence levels, we chose 7^=1.1 and /3t=1-58 
and box size=64 pixels. With these values, the super- 
sonic models were able to determine the Mach number 
regime with 98% accuracy, while the subsonic cases had 
an accuracy of 67% (see Figure fTT|). This is not surpris- 
ing however, since the moments are known to be a more 
robust measure of highly supersonic turbulence (Kowal, 
Lazarian & Beresnyak 2007). We plot the LOS Mach 
number map for model number 1 in the top panel of 
Figure [T2l with contours from the kurtosis moment map. 
The moments trace areas where the sonic Mach number 
is changing. 

We apply the moving box method to the SGPS data 
cuts using the same parameters for 77^, /3t and box size 
as were used for the simulations. For the SGPS data we 



obtained average skewness and kurtosis values for this 
moment map of 0.3 and 0.9 for skewness and kurtosis 
respectively. We over-plot the SGPS data with the kur- 
tosis moment map contours in the bottom panel of Fig- 
ure [121 While these values seem to point in the direction 
that the SGPS data is subsonic or transonic, we note 
that this method is less accurate for subsonic type tur- 
bulence. However, it does confirm that the gas in this 
patch of the sky is not statistically similar to what is 
expected for supersonic flows. 

4.2. Topology: The Genus Statistic 

The fllaments seen in |VP| show substantially different 
morphology when comparing maps of subsonic and su- 
personic turbulence, as was discussed in Section 3. Thus, 
a natural avenue of characterization would be to use 
topological measures in order to pick out different struc- 
tures. In this section we will investigate the utility of the 
genus statistics in order to characterize the topology of 
I VP I filaments. 

The genus statistics was developed to study the topol- 
ogy and deviations from Gaussianity of the universe and 
the distribution of galaxies in three dimensions (Gott 
et al. 1986; 1987). The use of genus statistics for the 
study of HI was first discussed in Lazarian (1999), and 
subsequent studies presented the genus curves for the 
SMC (Lazarian et al. 2002; Lazarian 2004, Chepurnov 
et al. 2008) and for MHD simulations (Kowal, Lazarian 
& Beresnyak 2007). 

Genus is a quantitative measure of topology. It can 
characterize both 2D and 3D distributions. Generally 
speaking, the genus is used to detect departures for Gaus- 
sianity. When dealing with the ISM one cannot expect 
deviations from symmetry to be small, especially in the 
presence of supersonic flows. In this case, genus can be 
used to characterize flows that are supersonic since these 
show large deviations from Gaussianity. 

The 2D genus can be represented as (Coles 1988; 
Melott et al. 1989): 

G = (# isolated high-density regions) 
— (# isolated low-density regions) 

where low- and high-density regions are selected with re- 
spect to a given contour threshold. For instance, a uni- 
form circle would have a genus of (one connected region 
of high density, i.e., an "island," and one connected re- 
gion of low density, while a ring (a donut, for example) 
would have a genus of -1 (one connected region of high 
density and two connected regions of low density) . Thus 
the genus can distinguish between "meatball" and "swiss 
cheese" topologies (Gott 1990). 

For a 2D image, the genus is simply a number which 
corresponds to a given threshold value. What is con- 
sidered to be high or low is dependent on the threshold 
value, which acts as a free parameter. As a result, for 
a given 2D image, the threshold value can be varied to 
construct a curve with the genus value on the y-axis and 
the threshold value on the x-axis. This is known as the 
genus curve. We show an example genus curve for |VP| 
in Figure 1131 This genus curve is for the subsonic sub- 
Aflvenic (model number 1) simulation. In the case of 
a density or column density field, the genus of subsonic 
turbulence is close to a Gaussian field. However the | VP| 
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Fig. 9. — Maps of |VP| for model 6 with smoothing FWHM=6 pixels (right) and without smoothing (left). The telescope resolution has 
a strong effect on the distribution of |VP| yet the morphology remains similar. 
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Fig. fO. — Moments of |VP| vs. smoothing for four different 
models. Error bars are created by taking the standard deviation 
between different time snapshots of the data. Smoothing FWFIM 
is given in pixels. 

distribution has longer tails and a more pronounced min- 
imum. The maximum and minimum points of the genus 
curve correspond to percolation of the distribution (see 
Colombi et al. 2000). The fact that the observed genus 
falls more slowly at large thresholds than the Gaussian 
distribution indicates that the islands are more discrete 
and pronounced than for the Gaussian distribution. 
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Fig. f 1. — The success rate of the moment map. A moment map 
is constructed with parameters 7^=1.1 and /3t=1.58 and box size 
(64 pixel), which were chosen with a genetic algorithm, ff a value 
in the moment map is above both 7^ and fix, then that pixel is 
considered supersonic. If the value is below either, then the pixel is 
considered subsonic. The moment map is compared with the actual 
LOS sonic Mach number map to determine whether the moment 
map alone is sufficient for recovering regimes of supersonic or sub- 
sonic on a pixel-by-pixel basis. The moment map is very successful 
(almost 100%) in determining if supersonic turbulence is present. 
The success rate drops to si60% in the case of diagnosing subsonic 
turbulence. This further confirms the higher order moments utility 
in the presence of supersonic flows. 



We expect that the sign of the genus curve at the mean 
intensity level of |VP| (i.e. where the curve crosses the 
X-axis) does describe the field topology. A positive genus 
will represent a clump-dominated field, while a nega- 
tive one will mean the domination of holes. However, 
it is more convenient to work with the zero of the genus 
curve (genus zero) vq, because it can be normalized to 
the field variance. In other words, consider a map in 
which the mean value of the map is subtracted off. vq 
represents the location where the genus curve crosses the 
X-axis, with the origin being the mean value of the field, 
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Fig. 12. — Top: LOS sonic Mach number map for model num- 
ber 1 with over-plotted contours of the kurtosis moment map with 
7t = 1.1 and /3t=1.58 and smoothing kornel= 64 pixels. The mo- 
ment map generally outlines areas where the sonic Mach number 
is changing, which is expected for a gradient quantity. Bottom: 
Map of the SGPS test region used in the paper with the kurtosis 
moment map contours over-plotted. 

SO that 0(1^0) — 0. In this case, for the intensity with the 
subtracted mean value, a negative vq corresponds to the 
clumpy topology, while a positive vq indicates the "swiss 
cheese" topology. 

We make plots of the genus zero vs. smoothing in Fig- 
ure [m A negative shift indicates a clumpy topology, 
while a positive shift indicates a hole topology. The er- 
ror bars are derived by estimating the variance of the 
genus distribution. For the error bar creation, we follow 
the method used by Chepurnov et al. (2008) in that we 
generated a set of images with randomly shifted phases 
of individual harmonics. This procedure causes the field 
to take Gaussian statistics, and therefore approximately 
the genus zero is at the origin. However, slight devia- 
tions from Gaussian allows us to effectively estimate its 
variance. The procedure is as follows: we take a fast 
Fourier transform (FFT) of the region being studied and 
assign the phase of each harmonic to a random variable 
uniformly distributed. After the inverse FFT, we calcu- 
late the respective Vq- After repeating this procedure ten 
times, we calculate the variance of the vq values. 

Figure [14] plots the genus zero vs. smoothing because 
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Fig. 13. — Genus curve of |VP| for model number 1 (pictured 
in the bottom right of Figure [!}• The genus curve of a Gaussian 
distribution is shown in dashed lines. The |VP| genus deviates 
strongly from the Gaussian model. Horizontal and vertical solid 
lines reference the origin. 
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Fig. 14. — Genus zero values of simulated | VP| . Negative values 
imply the topology is clumpy (meatball) while positive values imply 
the topology is hole dominated (swisscheese). The genus of | VP| is 
fairly insensitive to smoothing of maps of Q and U. The supersonic 
cases show a hole topology while the subsonic case shows a clump 
topology, even in the case of smoothing. Transonic cases show 
topology that is a mixture of clumpy and hole dominated. 

the telescope resolution is an important quantity to con- 
sider in any statistical analysis of observational data. As 
was demonstrated with the PDF moments in the previ- 
ous subsection, as maps of Q and U are smoothed, the 
distribution tends towards Gaussian. As such we include 
analysis of the genus zero with smoothing to see what 
effect the telescope resolution has on the topology. We 
include the genus zero vs. smoothing plots for maps of 
I VP I in the supersonic, transsonic and subsonic turbu- 
lence regimes. We see that over a range of smoothing 
values the topology between the three regimes is very 
different. The supersonic case shows a positive genus 
zero indicating that the topology is more hole like. The 
topology of the subsonic case is more clump like. Visual 
inspection of Figure [5] confirms this. In the case of the 
transonic turbulence, the topology is more neutral. In- 
terestingly, there is not a strong dependency when we 
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Fig. 15. — Genus curve for | VP| of the SGPS data set. The genus 
zero is at -0.07 which indicates a slightly clumpy topology. This 
is most similar to the transonic type genus curves {A4b ~ 2.0). 
Horizontal and vertical solid lines reference the origin. 

smooth the maps of Q and U. The topology remains 
roughly consistent over a range of smoothing. 

We show the genus curve for the SGPS |VP| data 
in Figure [151 The genus zero is at -0.07, which in- 
dicates a slightly clumpy topology. This value more 
closely matches the transsonic values accounting also for 
smoothing effects. 

5. DISCUSSION AND CONCLUSIONS 

Quantitative analysis (via PDF moments and genus) 
of the polarization gradient indicates that turbulence in 
the ionized ISM in these directions is in the range of 
subsonic to transonic. Our findings are supported by re- 
cent studies of Balmer-a emission measures, which have 
similarly found Mach number of 2 or less (see Hill et 
al. 2008). For studies of neutral gas, i.e. 21 cm HI 
gas, the warm phase of the SMC was shown to have 
properties of transonic gas (Burkhart et al. 2010). In 
the case of the cold component, several independent sta- 
tistical and direct measurements applied to the spectral 
lines and the column density indicate this component of 
the SMC is supersonic (Stanimirovic & Lazarian 2001, 
Burkhart et al. 2010). Furthermore, a recent study by 
Chepurnov & Lazarian (2010) added the Wisconsin H- 
a Mapper (WHAM) emission measure data to the "Big 
Power Law" of electron density fluctuations (see Arm- 
strong et al. 1995) provided further proof for a parsec 
to AU -5/3 power law, indicating that this phase of the 
ISM is not highly compressible. Thus it becomes appar- 
ent that different types of data, sampling different phases 
of the ISM, can compliment each other and provide more 
clues towards furthering our understanding of a turbu- 
lent ISM. 

Observational studies of ISM turbulence are extremely 
important. From the observational data one would like 
to obtain the characteristics of turbulence, e.g. intensity 
of its driving, importance of magnetic field in the dynam- 
ics of the media. These characteristics can be evaluated 
if we know sonic and Alfven Mach numbers, i.e. Ms and 
Ma- There have been several attempts to develop the 
techniques to get these numbers from observations using 
both column density data (see Kowal, Lazarian & Beres- 
nyak 2007 Burkhart et al. 2009) as well as spectroscopic 



data (Burkhart et al. 2011a). An example of success- 
ful application of such techniques to the Small Magel- 
lanic Cloud (SMC) HI data is provided in Burkhart et 
al (2010). 

This paper explores the utility of a new observational 
motivated technique, namely the gradient of the polariza- 
tion map, to determine Ms and Ma- The observational 
advantage of using gradients stems from the fact that 
these gradients are easily available from interferometric 
observations. Our study shows that |VP| is a very useful 
measure which allows one to study turbulence. We might 
also expect gradients be useful for turbulence studies in 
other types of data sets where shock morphology will be 
observed (i.e. column density, see Figures [2] and [3l) . 

We view this paper as a first taste of the utility of 
using polarization data (and their gradients) for studies 
of turbulence. There are many additional avenues that 
should be explored from this first step. For instance, 
what are the effects of multiple screens along the LOS? 
What are the effects of changing the assumption of the 
constant background polarization? What are the effects 
of other equations of state and the inclusion of partial 
ionization? 

One may wonder whether an additional way of study- 
ing turbulence is valuable, if we already have a few other 
ways to study turbulence, e.g. with column densities 
and position-position- velocity (PPV) spectroscopic data 
cubes. The answer is a sounding yes!. First of all, dealing 
with as complex media as the ISM, we would like to have 
as many independent measures as possible. Second, dif- 
ferent measures may be more sensitive to different phases 
of the ISM (see the list of the idealized phases and their 
magnetizations in Yan, Lazarian & Draine 2004). For 
example, the rotation measure, linear polarization and 
their gradients are biased towards ionized parts of ISM. 

The Alfven and sonic Mach number are not the only 
characteristics of the ISM turbulence. Spectra of turbu- 
lence, its intermittencies (see Kowal & Lazarian 2010) 
provide other measures which should be used to study 
interstellar turbulence in its complexity. 

5.1. Conclusions 

We created maps of the spatial gradient of the polariza- 
tion vector of isothermal MHD simulations and compared 
these with observations. We tested two statistical meth- 
ods on gradient polarization maps, namely the genus and 
higher order moments of the distribution to determine if 
these statistics were sensitive to the sonic Mach number. 
We found that: 

• Filamentary structure was created over a range of 
sonic Mach numbers, including cases where both 
shocks and subsonic turbulence were present. 

• Filaments showed different morphology for differ- 
ent regimes of sonic Mach number. 

• I VP I maps with high sonic Mach number showed 
filaments with a double jump profile that traced 
shocks, while subsonic cases showed filaments that 
were due to random fluctuations in the rotation 
measure along the LOS. Transonic simulations had 
both types of filaments present 
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• The moments of the |VP| distribution were higher 
for larger values of sonic Mach number but were 
also sensitive to the telescope resolution. 

• There is a strong Alfvenic dependency in the mo- 
ments of polarization angle and | VP| for sight-lines 
parallel to the ordered magnetic field. 

• The skewness and kurtosis moment maps of |VP| 
were successful at picking out subsonic or super- 
sonic pixels 67% and 99% of the time, respectively. 

• The genus of |VP| revealed a "hole" topology for 
supersonic cases and a "clump" topology for sub- 
sonic cases. Transsonic cases showed neutral topol- 
ogy. This trend is generally smoothing indepen- 
dent. 

• We applied the PDF moments and genus to the 
SGPS test region and found that this area was sta- 



tistically similar to models of subsonic to transsonic 
type turbulence. 
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APPENDIX 
THE GENETIC ALGORITHM 

We employ a genetic algorithm to search for the optimal three parameters (7t and /3t and box size) in the creation 
of the moment map, as discussed in Section 4. The literature on the genetic algorithm is vast, so we refer the curious 
reader to Eiben & Smith (2003) and references therein. Unlike brute force methods, the genetic algorithm relies on 
principles of biological evolution, such as reproduction and natural selection, to determine the optimal parameters. 

A typical genetic algorithm requires a fitness function to evaluate the solution domain. In our case the fitness metric 
is how well the local skewness and kurtosis arc able to determine if the local sonic Mach number is either in the 
supersonic or subsonic regime. In our case, local means on a pixel-by-pixel basis. The main advantage of us choosing 
a genetic algorithm over a brute force method is the cost in computation time. Rather then looping over all possible 
combinations of parameters, we only seek a convergence to a range of the best fit outcomes. We ran the algorithm 
several separate times on multiple simulations with varying sonic and Alfvenic Mach numbers in order to assure that 
the convergence obtained was repeatable. 

The algorithm we employ is a simple implementations of a genetic algorithm and works as follows. We start with a 
population of box sizes and threshold values. The algorithm initially picks a random subset of ten candidates from this 
larger population pool; that is, it picks ten combinations of 7t, I3t, and box size. Then the moment map is created by 
calculating the moments in the kernel, then moving the kernel pixel-by-pixel and repeating the moment calculation. 

We take this map and compare it to the LOS sonic Mach number map (smoothed by the same degree at the moment 
map) pixel-by-pixel. If a pixel in the moment map is above "fx and /3t then it is considered supersonic, else it is 
subsonic. If this matches with what is seen in the LOS sonic Mach number map then we give this pixel a value of 1, 
if it does not match then we give this pixel a value of 0. We calculate the percentage of successful pixels to determine 
the fitness of that particular choice of parameters. 

Once we have determined the percent of success for all ten of these initial candidates, we "clone" the most fit of 
these and then "breed" the candidates to fill out the other nine models of the next generation. This keeps the total 
population size constant. Our breeding is done by averaging the parameters from a random selection of these ten 
candidates (the parent population). The most successful candidate from the parent population and the averaged 
candidates now become the "children population." The algorithm is repeated from the point of making the moment 
map. This process continues for 10 iterations of parent and child. Because our algorithm uses small population 
sizes, after a fairly small number of generations the variety in free parameters will shrink significantly, resulting in an 
"inbred" population. Therefore every 10th generation, instead of breeding the next generation, we clone the fittest 
model and then fill out the rest of the subsequent generation with new models, randomly chosen in the same way as 
the models used in the first generation of the algorithm. This process ensures that the genetic algorithm always has a 
large region of parameter space to explore. 

We ran the genetic algorithm multiple times to ensure the solution of best fit was converging to roughly the same 
parameters. While there were many parameter sets in the 90%-|- success range in the supersonic case, we choose 
7t=1.1 and /3t=1-58 and box size=64 pixels, as it was the best match we found with the subsonic cases (67%) 
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